#include "matrix.h"

/* Changes the internal memory organization. This should only be touched for performance 
   reasons, if at all. That's why this is not exposed in util.h */
int row2colMajor(matrix *A) {return _ENOTIMPL;}
int col2rowMajor(matrix *A) {return _ENOTIMPL;}

int init_mat(matrix *A, int row, int col) {

    double *d = (double *)malloc(row*col*sizeof(double));

    if (d == NULL) 
        return _ENOMEM;

    A->row = row;
    A->col = col;
    A->major = ROWMAJOR;    /* Let's make 'em row major by default */
   
    A->data = d;

    return _ESUCCESS;
}

void dest_mat(matrix *A) {
    free(A->data);
}

int add_mat(matrix *A, matrix *B, matrix *res) {

    int rv;
    int row = A->row;
    int col = A->col;

    /* Check both matrices have the same size */
    if (!(A->row == B->row && A->col == B->col)) 
        return _EBADSIZE;

    /* initialize the result matrix */
    rv = init_mat(res, row, col);

    if (rv != _ESUCCESS)
        return rv;


    for (int i=0; i<row; i++) {
        for (int j=0; j<col; j++) 
        /* Do the actual addition */            
        MAT_ELEM(res,i,j) = MAT_ELEM(A,i,j) + MAT_ELEM(B,i,j);
    }

    return _ESUCCESS;
}

int mult_scalar(matrix *A, double t, matrix *res) {

    int rv;
    int row = A->row;
    int col = A->col;

    /* initialize the result matrix */
    rv = init_mat(res, row, col);

    if (rv != _ESUCCESS)
        return rv;

    for (int i=0; i<row; i++) {
        for (int j=0; j<col; j++) 
        /* Do the actual multiplication */            
        MAT_ELEM(res,i,j) = MAT_ELEM(A,i,j) * t;  
    }

    return _ESUCCESS;
}

int mult_mat(matrix *A, matrix *B, matrix *res) {

    int rv;
    int m = A->row;
    int n = A->col;     // n = B->row too (if sane input)
    int p = B->col;

    if (A->col != B->row) 
        return _EBADSIZE;

    /* initialize the result matrix */
    rv = init_mat(res, m, p);

    if (rv != _ESUCCESS)
        return rv;


    /* Calculate the (i,j)th element of res */
    for (int i=0; i<m; i++) {
        for (int j=0; j<p; j++) {

            double temp = 0.0;
            for (int k=0; k<n; k++) {
                temp += MAT_ELEM(A,i,k) * MAT_ELEM(B,k,j);
            }
            MAT_ELEM(res,i,j) = temp;
        }
    }

    return _ESUCCESS;
}

int trans_mat(matrix *A, matrix *res) {

    int rv;
    int row = A->row;
    int col = A->col;

    rv = init_mat(res, col, row);

    if (rv != _ESUCCESS)
        return rv;

    for (int i=0; i<row; i++) {
        for (int j=0; j<col; j++) {
            MAT_ELEM(res,j,i) = MAT_ELEM(A,i,j);
        }
    }

    return _ESUCCESS;
}
